clc;clear;close all
%% 读取数据预处理
data=readtable('负荷天气数据.csv','range','B2:B10001');
data=table2array(data);
[x,y]=data_process(data,30);%前30个时刻 预测下一个时刻
% 归一化
[xs,mappingx]=mapminmax(x',0,1);x=xs';
[ys,mappingy]=mapminmax(y',0,1);y=ys';
% 划分数据
n=size(x,1);
m=round(n*0.7);%前70%训练，对最后30%进行预测
XTrain=x(1:m,:)';
XTest=x(m+1:end,:)';
YTrain=y(1:m,:)';
YTest=y(m+1:end,:)';
%% 设计网络
numFeatures = size(XTrain,1);% 输入特征数
numResponses = 1;% 输出响应数
numHiddenUnits = 200;% 隐含层数目
%  GRU设计
layers1 = [
    sequenceInputLayer(numFeatures,"Name","sequence")
    gruLayer(numHiddenUnits,"Name","gru")
    fullyConnectedLayer(numResponses,"Name","fc")
    regressionLayer("Name","regressionoutput")];
% LSTM设计
layers2 = [ ...
    sequenceInputLayer(numFeatures)
    lstmLayer(numHiddenUnits)
    fullyConnectedLayer(numResponses)
    regressionLayer];
% 训练优化设计
options = trainingOptions('adam', ...% adam形式梯度下降优化
    'MaxEpochs',200, ...% 最大迭代次数
    'GradientThreshold',1, ...
    'InitialLearnRate',0.001, ...% 初始学习率
    'LearnRateSchedule','piecewise', ...
    'LearnRateDropPeriod',125, ...
    'LearnRateDropFactor',0.2, ...
    'Verbose',0);
% 训练两种方案网络
net1 = trainNetwork(XTrain,YTrain,layers1,options);
net2 = trainNetwork(XTrain,YTrain,layers2,options);
numTimeStepsTest = size(XTest,2);
% 对网络1结果分析
for i = 1:numTimeStepsTest
    [net1,YPred(:,i)] = predictAndUpdateState(net1,XTest(:,i),'ExecutionEnvironment','cpu');
end
% 反归一化
predict_value=mapminmax('reverse',YPred,mappingy);
true_value=mapminmax('reverse',YTest,mappingy);
disp('结果分析')
rmse=sqrt(mean((true_value-predict_value).^2));
disp(['根均方差(RMSE)：',num2str(rmse)])
mae=mean(abs(true_value-predict_value));
disp(['平均绝对误差（MAE）：',num2str(mae)])
mape=mean(abs(true_value-predict_value)/true_value);
disp(['平均相对百分误差（MAPE）：',num2str(mape*100),'%'])
fprintf('\n')
figure
plot(true_value,'-*','linewidth',1)
hold on
plot(predict_value,'-s','linewidth',1)
legend('实际值','预测值')
grid on
title('GRU')
